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A quantum Monte Carlo algorithm for the transverse Ising model with arbitrary short- or long- 
range interactions is presented. The algorithm is based on sampling the diagonal matrix elements 
of the power series expansion of the density matrix (stochastic series expansion), and avoids the 
interaction summations necessary in conventional methods. In the case of long-range interactions, 
the scaling of the computation time with the system size N is therefore reduced from N 2 to N In (N). 
The method is tested on a one-dimensional ferromagnet in a transverse field, with interactions 
decaying as 1/r 2 . 
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I. INTRODUCTION 

Monte Carlo studies of classical and quantum many- 
body systems with long-range interactions are limited by 
time-consuming summations over the interacting parti- 
cle pairs, the number of which grows quadratically with 
the system size. Many important problems in both basic 
and applied science can be mapped onto long-range in- 
teracting spin models, and hence it would be desirable to 
develop more efficient numerical techniques for tackling 
them. For classical Ising models, considerable progress 
has indeed been made on algorithms scaling almost lin- 
early with the system size . In the context of simulated 
annealing 0, where the ground state of a classical sys- 
tem (typically with complicated interactions) is obtained 
through a simulation where the temperature is slowly 
lowered to zero, it has been suggested |3j that a more 
rapid convergence could be achieved by using a quantum 
model, e.g., the Ising model in a transverse (spin flipping) 
field. Even in an imaginary-time path-integral formula- 
tion, the quantum fluctuations can, at least in some cases 
0, relax the system towards its classical ground state 
more rapidly than thermal fluctuations. This is a strong 
motivation for developing more efficient simulation meth- 
ods for quantum Ising models. Another important reason 
is the continued prominence of the transverse Ising model 
in the theory of magnetism, particularly in the context of 
quantum phase transitions [a, IE • Whereas transverse 
Ising models with short-range interactions have recently 
been actively studied using quantum Monte Carlo meth- 
ods 0, , numerical work on long-range models has so 
far been limited to special cases j£f ■ In some of the best 
experimental realizations of the transverse Ising model 
the interactions are in fact long-ranged 0- 

Here a stochastic series expansion (SSE) [13 algorithm 
for transverse Ising models with long-range interactions 
is introduced in which the direct summation over the in- 
teracting spins is avoided. The computation time scales 
with the system size N as N In (N) times the spatial in- 
tegral of the absolute value of the interaction [which nor- 
mally converges as N — > oo, or diverges only as ln(.ZV)]. 



Both local and cluster-type updates are developed for 
the transverse Ising model with arbitrary interactions. 
The cluster update is a generalization of the classical 
Swendsen-Wang cluster method to the transverse 
Ising model, and shares some features with a scheme pre- 
viously used within the continuous-time world-line algo- 
rithm 0. The way to treat the long-range interactions 
generalizes the scheme developed for the classical Ising 
model by Luijten and Blote 0,0. The integration of 
these features in the SSE formalism should open new op- 
portunities for detailed numerical studies of a wide range 
of important models. The algorithm is here tested on a 
ferromagnetic chain with interactions decaying as 1/r 2 , 
for which results in the classical limit are available for 
comparison [13 . fpiL fl5L 1161 1X711 . 

In Sec. II the application of the SSE method to the 
transverse Ising model is described in detail. Local up- 
dates as well as classical and quantum-cluster updates 
are discussed. Results for the model with 1/r 2 interac- 
tions are presented in Sec. III. Sec. IV concludes with a 
brief discussion. 



II. STOCHASTIC SERIES EXPANSION 

The SSE method jn| is an efficient alternative to 
worldline quantum Monte Carlo 0|. It is based on a 
generalization of the power-series scheme for the Heisen- 
berg ferromagnet that was developed by Handscomb in 
the early 1960's [l^. Handscomb's method was later ex- 
tended to some other models [2(J, but the requirement 
of analytically calculable traces of the terms of the ex- 
pansion inhibited further progress. In the SSE method, a 
basis is instead chosen, and the traces are also evaluated 
stochastically, in combination with the sampling of the 
operator products in the series expansion of exp(— /3H). 
This starting point for quantum Monte Carlo is as gen- 
erally applicable as the worldline (imaginary-time path- 
integral) approach. Recently, loop-type cluster updates 
[2l| have been developed and generalized for efficient SSE 
simulations of a wide range of models [22ll23l|. However, 
since the loop updates rely heavily on the presence of off- 
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diagonal pair (or multi-particle) interactions, they can- 
not be directly adapted to the transverse Ising model in 
the standard basis where the Ising term is diagonal. In 
the basis where the field is diagonal, loop updates can be 
easily implemented [2^, but then sign problems [24| 
appear when the interaction is frustrated. Here the SSE 
method is applied to an arbitrary transverse Ising model, 
i.e., with no limitations on the sign and range of the spin- 
spin interaction. Several types of local and cluster-type 
updates will be described. 

A. Configuration space 

Consider the general Hamiltonian for the Ising model 
in a transverse field of strength h, 

where Oi is a Pauli spin matrix (<r| = ±1) and is 
the strength of the interaction between spins i and j, 
which can be random or uniform and of any sign. The 
dimensionality is arbitrary. Define the operators 

Ho,o = 1, (2) 

H it o = h{a++a~), z>0, (3) 

H hl = h, z>0, (4) 

Hij = \Jij\ - Jij-ofoJ, i,j > 0, i ^ j. (5) 

Up to a constant, the Hamiltonian can be written as 

N N 

II ( 6 ) 

»=1 3=0 

The constants Hij are introduced for purposes that will 
become clear below. Note that 7?o,o is not included as a 
term in the Hamiltonian © but will be important in the 
simulation scheme. 

In the SSE approach |Toj to finite-temperature 
quantum Monte Carlo, the partition function Z = 
Tr{exp(— /3H)} is written as a power-series expansion, 
with the trace expressed as a sum over diagonal matrix 
elements in a suitably chosen basis. Using |JB} then gives 

oo an n 

3 = £E£^< a in*««wwi«>. ( ? ) 

a n=0 S n ' 1=1 

where S n denotes a sequence of n operator-index pairs 
(hereafter referred to as operators): 

S n = [i(l),j(l)},...,[i(n),j(n)}, (8) 

with e {1, . . . , N} and £ {0, ... , TV}. The stan- 
dard basis {\a}} = {|of, ■ ■ ■ ,0^)1 is used. 

Because of the constants added to Hij in the 
eigenvalues of these operators are 2|Jy| and 0. All non- 
zero terms in (JJJ are therefore positive and can be used as 



relative probabilities in an importance sampling scheme. 
A term is specified by a state \a) and an operator se- 
quence S n - One can show that the total internal ene rgy 
(including the constants added to H) is given by pi Eg 
E = —(n)/(3. Hence, the size of the operator sequence to 
be stored in computer memory scales as /3NIn(J), where 

N N 

WH-tfEEi-u ( 9 ) 

i=l i=i 

which converges or grows much slower than N for most 
cases of interest. 

In order to construct an efficient sampling scheme, it 
is useful to cut the expansion Q at some power n = 
L, sufficiently high for the remaining truncation error 
to be exponentially small and completely negligible [L 
clearly has to be ~ /3NIn(J)]. One can then obtain an 
expansion for which the length of the operator sequence 
is constant, by considering random insertions of L — n 
unit operators Hq q in the product in (J7J. Adjusting for 
the ( ) possible insertions gives 

1 L 

^ = ^EE W - ™) ! ("I II H m,m l«> . (io) 

a S L 1=1 

where — [0,0] is now also an allowed operator 

in the sequence Sl, and n denotes the number of non- 
[0, 0] operators. Note again that Hq o is not part of the 
Hamiltonian, but is introduced only for the purpose of 
constructing a computationally simpler updating scheme 
where the operator list has a fixed length. 

It is useful to define states \a(p)) = \crf(p), . . . , <J z N {p)) 
obtained by propagating \a) — |ct(0)) by the first p oper- 
ators in Sl] 

v 

|a(p)> = r Y[ H m,JQ) l a )' (n) 

where r is a normalization factor. A non- vanishing ma- 
trix element in l|10(l then corresponds to the periodicity 
condition \a{L)) = |a(0)), which requires that for each 
site i there is an even number (or zero) of spin flipping op- 
erators [i, 0] in Sl- Definition (JJJ implies that the Ising 
operators may act only on states with of = cr| if 
Jij < (ferromagnetic), or erf = —a* if Jy > (antifer- 
romagnetic). There are no other constraints. 

An SSE configuration is illustrated in Fig. The ver- 
tical direction in this representation will be referred to 
as the SSE propagation direction. It can be related to 
the imaginary-time direction in standard path integral 
representations |2fl| . Note that this full configuration, 
including all the states \a(p)) explicitly, does not have 
to be stored in the simulation. A single state and the 
operator sequence suffice for reproducing all the states, 
and such a representation is used in some stages of the 
simulation. For some updates it is convenient to generate 
other representations, as will be discussed below. 
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FIG. 1: An SSE configuration for an 8-site one-dimensional 
system. Here the truncation L — 49, and the expansion or- 
der of the term (i.e., the number of Hamiltonian operators 
present) n = 40. The solid and open circles represent the 
spins c|(p) = ±1, with the propagation index p = 0, . . . ,L 
corresponding to the different 8-spin rows. The thick and thin 
short horizontal bars represent spin-flip operators Hi : o and 
constants Hi t i, respectively. The longer lines represent Ising 
operators Hij (i ^ j) acting on the spins at the line-ends. 



B. Local updates 

The sampling of Eq. (|10|) can be carried out using sim- 
ple operator substitutions of the types 



[0,0], 
[i, 



[i,0] Pl [i,0] P2 , 



(12) 
(13) 



where the subscript p indicates the position (p = 
1, . . . , L) of the operator in the sequence Sl- The power 
n is changed by ±1 in the diagonal update (I12fl and is 
unchanged in the off- diagonal update (|13l) . In the di- 
agonal update the Ising terms [i,j] and the constants 
are sampled. The constants are used in the off- 
diagonal update as a means of achieving easy insertions 
and removals of two spin-flipping operators [z, 0]. With 
the value h chosen for the constant in Q, the operator 



replacements do not change the weight of the SSE con- 
figuration. However, the off-diagonal update also leads 
to spin flips in the propagated states between p\ and 

P2] of Oi), ■ ■ .,0?(p 2 - 1) -> -of(Pi), ■ ■ ■ , -0?(p 2 - !)• 
[pi > Pi also has to be considered, leading to flipped 
of (pi), . . . , of (L - 1) of (0), . . . , of (p 2 - 1)], which is al- 
lowed if (and only if) no Ising operators acting on site 
i are present in Sl between positions p\ and pi. Note 
that this constraint is completely local, regardless of the 
range of the interaction, and that the update requires no 
knowledge of the spin state. This is the reason for the ad- 
vantage of this simulation scheme over worldline methods 
Il8| - where calculating the acceptance probability for 
every update requires a summation over all the spins in- 
teracting with those flipped. Here an allowed off-diagonal 
update l|13|) leaves the weight unchanged and can be car- 
ried out with probability one. 

If h ^ 0, the above updates of the operator sequence 
suffice for achieving ergodicity. If there are no Ising op- 
erators acting on a site i, af (0), . . . , of (L — 1) can also 
be flipped without changes in Sl- This update in princi- 
ple makes simulations using the present scheme possible 
also for h = 0, but in practice unconstrained spins occur 
frequently only at high temperatures, when (n) is small. 
Other types of "classical" spin flips — flips of clusters — 
are also possible, and will be discussed in Sec. II C. 

The simulation can be started with a random state 
a(0)) and a sequence Sl containing only [0, 0] operators. 
The truncation L can be chosen arbitrarily (small); it is 
adjusted during the equilibration part of the simulation, 
e.g., by requiring L > (4/3)n after each update. This en- 
sures than n never reaches L during the reminder of the 
simulation, and hence that there will be no detectable 
systematic errors arising from the truncation of the ex- 
pansion |l0j| . In the beginning of an updating cycle, the 
operator sequence Sl and the state |a(0)) is stored. 

The diagonal update (|12|) is attempted successively for 
all p = 1, . . . , L. In the course of this process, the spin 
state is propagated by flipping spins of as off-diagonal 
operators [i, 0] are encountered in Sl, so that the states 
\a(p)) are generated successively. For an — > [0,0] 
update, i.e., removing a Hamiltonian operator, there 
are no constraints and the update should always be ac- 
cepted with some non-zero probability. In the case of 
[0,0] — > i.e., inserting an operator from the Hamil- 
tonian, there are constraints, and the update may not be 
allowed for all However, initially the indices i,j are 
left undetermined and it is assumed that any would 
be allowed. Under this assumption, the acceptance prob- 
abilities for the diagonal update are given by 

P(Nh 



P([0,0] 
P([i,j] 



•M) = 

[0,0]) = 



2£«|J«I 



L 



— n 
n + 1 



f3{Nh - 



2E«I^-I 



(14) 
(15) 



where Y]fj does not include i = j and P > 1 should be 
interpreted as probability one, as usual. These probabil- 
ities are simply obtained from the ratio of the new and 
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old prefactor in l|10|) when n — ► n ± 1 ; 



±1 [L-(n±l)]! 



(L-n)\ 



(16) 



and the ratio between the matrix element 1 of the [0, 0] 
operator and the sum Nh + 2 I Jij I °f the non-zero 
matrix elements of all operators. Staying with the 
assumption that any is allowed in the update [0, 0] — > 
the relative probability of an operator with the first 
index i is P(i) = ^2jMij, where My is the non-zero 
matrix element corresponding to H%j (i.e., h for i = j and 
and 2 1 Jy| else). The normalized cumulative probabilities 
P c (k — 1, . . . , N) are stored in a pre-generated table; 



P c (k) 



(17) 



In order to select the first index i of the operator to 
be inserted, a random number < R < 1 is generated. 
The table P c is searched (using, e.g., a simple binary 
search) for the smallest k for which P{k) > R; the first 
index of the operator [i,j] is then i — k. The second 
index can be chosen in a completely analogous way, with 
the relative probability for j, given i, being My. For a 
random system with long-range interactions, a pregen- 
erated table with iV 2 elements is hence needed for stor- 
ing all the cumulative probabilities for the second index. 
For non-random interactions in a translationally invari- 
ant system, the first index can be selected at random with 
equal probabilities without searching a table, and the size 
of the second table is reduced to N. For a short-range 
or truncated interaction the table size is smaller, corre- 
sponding to the number of spins within the range of the 
interaction; clearly, the whole selection process should 
then be reduced to a single step for obtaining both i and 
j (e.g., selecting one out of a total number ~ N of op- 
erators and reading the corresponding i,j from a table). 
The two-step procedure is advantageous for non-random 
long-range interactions, where it allows for the reduction 
of the size of the probability table from iV 2 to N. For 
random models, the storage requirement is always TV 2 , 
and it may then again be better to combine the first and 
second index searches, using a single size-iV 2 table for 
all the cumulative probabilities of For short-range 

random interactions the size of the table is N times the 
number of spins within the interaction range. 

The operator generated as above may or may not 
be allowed in the current spin configuration \a(p)}. If 
crf(p) and <Jj(j>) indeed are in an allowed state, is 
inserted at position p. Otherwise, the process for generat- 
ing [i, j] has to be repeated, until an allowed operator has 
been generated. The reject-and-repeat step leads to the 
correct probabilities for selecting among all the allowed 
diagonal operators Typically, an allowed operator 

is generated very quickly, since the interactions favor the 
allowed spin alignment. Note that the constants [i, i] are 
always allowed (for h > 0), so there is no risk of the 
search never terminating. 



The off-diagonal update l|13(l can be efficiently carried 
out if Sl is first partitioned into separate subsequences 
for each site i. Subsequence i contains only spin- flipping 
operators [i,0] and constants Their positions in Sl 
are also stored, to be used for recombining the subse- 
quences after the update. The constraints on modifica- 
tions at site i imposed by Ising operators or 
(for any j) can be stored as flags indicating the pres- 
ence of one or several of these operators between neigh- 
boring subsequence operators. Updating a subsequence 
amounts to selecting two non-constrained neighboring 
operators at random from the subsequence, and carry- 
ing out the substitution (| 1 31) if the two operators are 
identical. If they are different, they can be permuted. A 
number proportional to the subsequence length of such 
pair updates are carried out for each subsequence, after 
which they are recombined into a new Sl- 

The diagonal update H12|) at all positions in Sl require 
~ Lin (AT) ~ (3 N In (TV) I n{J) operations, where the fac- 
tor In (N) is the scaling of the average number of oper- 
ations needed to search the cumulative probability ta- 
ble(s) in the case of long-range interactions. Partitioning 
Sl into subsequences and updating all of them accord- 
ing to l|13f) requires on the order of L operations. Hence, 
the number of operations for a full updating cycle of the 
degrees of freedom of the system (one Monte Carlo step) 
scales as (3N In (N)Ijy( J). This should be compared to 
the (iN 2 scaling in worldline methods 0, 0] , where one 
power of iV is due to the summation required to calculate 
the weight change when flipping a spin interacting with 
N other spins. Here this summation has been circum- 
vented by writing the interactions in the SSE formalism 
as fluctuating constraints that are purely local. 



C. Classical cluster update 

In the Swendsen-Wang cluster algorithm for the 
classical Ising model, i.e., with h — and a uni- 
form nearest-neighbor interaction of strength J, auxiliary 
bond variables fry are introduced in order to construct 
clusters of spins that can be flipped independently of each 
other. Given a spin configuration, and with initially all 
bond variables &y = 0, for every interacting spin pair for 
which (Tj(jj = —J/\J\ (i.e., the orientation energetically 
favored) the bond variable is set, bij = 1, with probabil- 
ity P — 1 — e -2 '* 7 !' 3 . When all bonds have been visited, 
clusters of spins connected by bij = 1 bonds are formed, 
and each of these clusters is flipped with probability 1/2. 
Single spins not connected to any bij = 1 bond are single- 
spin clusters. After the clusters have been flipped, all the 
bond variables are again set to zero and the process is re- 
peated. This scheme can in fact be constructed using the 
SSE formalism, as an alternative to the Fortuin-Kastelin 
mapping |26| , on which the Swendsen-Wang algorithm is 
based. 

The relation to the Swendsen-Wang algorithm is shown 
as follows, by applying the SSE method to the classi- 
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cal Ising model, now again considering a general form 
of the interaction and with the bond-operator ify = 
| Jij | — Jij<Ji<Jj as in Eq. Since all operators com- 
mute, the operator e~P H can be written as a product of 
operators e@ Hi i = 1 + [3 Hij + . . .. The uniqueness of 
the power-series expansion then implies that in the SSE, 
where e~@ H is expanded directly, the probability of hav- 
ing one or more operators on a bond i, j when <7j<Xj = 
—Jij/\Jij\ is 1 — e _2 l Jij l' 3 , i.e., exactly the probability of 
having the bond variable bij = 1 in the Swendsen-Wang 
scheme. In a configuration a%aj — Jij/\Jij\ there can be 
no operators on the bond in the SSE, and the Swendsen- 
Wang = 1 probability is also zero per construction. 
One can hence make the connection that one or more 
operators acting on a spin pair in the SSE scheme corre- 
sponds to a filled bond (b^ = 1) in the Swendsen-Wang 
algorithm. The definition of a cluster is then exactly the 
same in the two algorithms. Clearly, such a cluster in the 
SSE can also always be flipped, since the Ising operators 
only impose constraints on the relative orientations of 
connected spins, which is maintained when the cluster is 
flipped. Since the weight does not change, the flip should 
be done with probability 1/2. The scheme is hence iden- 
tical to the Swendsen-Wang algorithm, except that the 
filled bonds bij = 1 in SSE are generated in a different 
way, using the diagonal update (fT^I . Note that for a 
classical model, all the propagated SSE states are 
identical, i.e., of (p) = of (0) for all p = 0, ...,L — 1, and 
hence no state propagations have to be considered as the 
diagonal update is carried out. 

It is interesting to note that the SSE scheme for the 
classical Ising model should in fact be more efficient than 
the standard Swendsen-Wang algorithm at high temper- 
atures. This is because the number of operators in the 
SSE operator list scales as E(T)/T, where E(T) is the 
total energy at temperature T (E ~ N) and for large 
T the construction of the clusters based on the operator 
list should then be faster than visiting all the bonds, as 
is done in the Swendesn-Wang algorithm. However, in 
practice the interesting physics occurs when the number 
of SSE operators per interacting spin pair is of the or- 
der of one or larger, and then there are no advantages of 
the SSE classical cluster algorithm relative to Swendesen- 
Wang. 

The classical SSE cluster update can also be used in 
the presence of the transverse field (h > 0). The clus- 
ters are defined in terms of bonds signifying the presence 
of one or more Ising operator, as above, without regard 
for the single-spin flipping operators Hi o and constants 
Hi i. These operators can be neglected because when a 
cluster is flipped, all spins of belonging to the cluster 
are implicitly flipped in all propagated states i.e., 
a i(p) ~ * ~ a i(p) for all p = 0, . . . , L — 1 (this is the rea- 
son for the term "classical cluster" even when h > 0) 
and hence all operations with the single-spin operators 
remain valid and produce the same factors in the weight 
before and after the cluster flips. Note again that only 
the first state, i.e., of (0), i = 1, . . . ,N, has to be stored 
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FIG. 2: Upper panel: Interaction bonds in a configuration for 
a 2D system with long-range interactions. Lower panel: The 
clusters constructed from the bonds. Sites with equal symbols 
belong to the same cluster. Dots indicate spins not acted on 
by any Ising operator and constitute single-spin clusters. 



when constructing the classical clusters. 

In the case of long-range interactions, a cluster can 
consist of several intertwined pieces on the lattice, as il- 
lustrated for a two-dimensional case in Fig. El Regardless 
of the range of the interaction, the construction of the 
clusters, given an SSE operator list, can be easily carried 
out using a number of operations scaling as the number 
of operators in the list. 

Since the classical SSE cluster update is equivalent to 
the Swendsen-Wang algorithm in the classical limit and 
only takes the Ising terms into account also in the quan- 
tum case, it cannot be expected to be efficient much be- 
yond the classical limit h — 0. For a non-random system 
that undergoes a phase transition at T c (0) when h = 0, 
the critical temperature is reduced by the transverse field; 
T c {h) < T c (0). Hence, the classical clusters will percolate 
for T > T c and this update will not be efficient close to 
T c . The primary reason to introduce the classical cluster 
update here was to demonstrate the relationship between 
SSE and the Swendsen-Wang algorithm. In the case of 
long-range interactions, the scheme becomes very similar 
to the Luijten-Blote algorithm [l|, again just differing in 
the way the bonds are generated. 
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D. Quantum-cluster update 

The purpose of the quantum-cluster update is to effect 
flips of spins of (p) only in a limited number of propagated 
states p, in different states for different sites i. In other 
words, these clusters will be finite and irregularly-shaped 
both in the space and SSE propagation (imaginary time) 
direction. In the process, operator substitutions Hij <-> 
Hifi (constant to spin-flip, and vice versa) will also be 
accomplished. This update hence replaces the local off- 
diagonal update l|13(l . 

To discuss the quantum-cluster u pda te, it is useful to 
introduce the notion of vertices [22L |23|| . Looking at the 
graphical representation of a configuration in Fig. ^ it 
can be noted that the vertical "lines" of same spins be- 
tween two operators acting on a given site constitute 
redundant information. The full configuration can be 
represented by a list of positions (on the lattice) of the 
operators, and the spin states (on one or two sites for 
the model considered here) before and after the opera- 
tors act. These relevant spins are called legs of the 2-spin 
vertices (corresponding to constant and spin-flip oper- 
ators) or 4-spin vertices (corresponding to Ising bond- 
operators). All possible vertices for the transverse Ising 
model are shown in Fig. [3] Note that only those Ising ver- 
tices that are compatible with the sign of the interaction 
between a given pair of spins are allowed for those spins; 
again, this is due to the choice of constant in the bond- 
operator JSJ. In the computer, the vertices are linked 
to each other by pointers, so that from a given vertex- 
leg one can reach the next or previous vertex that has 
a leg on the same site (i.e., there are links that replace 
the segments of vertical lines of same spins in Fig.^J. A 
detailed discussion of the practical implementation of a 
linked vertex list has been given in Ref. I23L 

To construct and flip a quantum-cluster, one of the 
legs of one of the n vertices is picked at random, and the 
corresponding spin is flipped. Depending on the type of 
the vertex, different actions are taken, examples of which 
are given in Fig. 0] The arrow pointing into the vertex 
indicates the entrance leg. In the case of an Ising vertex, 
all the four spins are flipped and the cluster building 
process branches out from all the legs, as indicated by the 
arrows pointing out from the vertex. Using the pointers 
of the linked vertex list, the arrows point to legs of other 
vertices; these become new entrance legs which are put 
on a stack and subsequently processed one-by one. If 
the entrance leg is on a constant or spin-flip vertex, only 
the entrance spin is flipped. The vertex type then also 
changes, in terms of operators from Hi^o to and 
vice versa. In these cases there is no branching-out and 
no new legs are put on the stack, i.e., this particular 
branch of the cluster terminates. If a link points to a 
spin that has already been flipped (i.e., two arrows point 
toward each other), that leg should not be used again as 
an entrance and is hence not put on the stack. Therefore, 
each vertex-leg can be visited at most once (each spin can 
be flipped at most once) and the cluster is completed 
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FIG. 3: All the possible 4-leg and 2-leg vertices, (a) Ferro- 
magnetic Ising vertices, (b) antiferromagnetic Ising vertices, 
(d) constant vertices, and (d) spin-flip vertices. 
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FIG. 4: Examples of vertex processes: (a) reversal of a ferro- 
magnetic Ising vertex, (b) constant to spin-flip, and (c) spin- 
flip to constant. 



when there are no more entrance-legs on the stack. The 
reason that the cluster can always be flipped is again 
that the SSE weight is not affected; the matrix element 
of the Ising bond-operator is not affected when both spins 
are flipped (in the absence of an external field in the in- 
direction, which would necessitate a modified approach), 
and the matrix elements for the constant and spin-flip 
operators are both equal to h. 

The construction of a single cluster, which is flipped 
with probability one, is a quantum-mechanical analogue 
of the classical Wolff algorithm |53 ; in the absence of the 
transverse field the clusters are identical to those of the 
Wolff algorithm. Note, however, that there is a difference 
when constructing more than one cluster: The number of 
operators in the SSE operator list, and their positions on 
the lattice, do not change in the quantum-cluster update. 
The clusters are therefore completely deterministic once 
the operator list is given. Hence, when constructing sev- 
eral clusters using the same SSE operator list, it is quite 
likely that the same cluster is constructed and flipped 
multiple times. This is clearly not efficient. However, one 
can also construct all clusters, as in the Swendsen-Wang 
scheme, and only flip them with probability 1/2. This is 
done by always starting a new cluster from a vertex-leg 
which has not yet been visited. Every vertex-leg belongs 
uniquely to one cluster, and clearly the number of opera- 
tions required to complete this updates then scales as L, 
i.e., typically as f3N. 

A natural definition of a Monte Carlo step including 
the quantum-cluster update is a full sweep of diagonal 
updates, followed by the construction of the linked list of 
vertices, in which all clusters are constructed and flipped 
with probability 1/2. After that, the updated vertex list 
is mapped back into a state \a(0)) and an operator se- 
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quence Sl- Free spins, i.e., those that are not acted on 
by any operators, can again be considered as single-spin 
clusters and should also be flipped with probability 1/2. 
No local off-diagonal updates l|13l) are needed. 

Since the quantum-cluster update explicitly includes 
the quantum mechanical features of the configurations 
(i.e., the presence of spin- flip operators), it can be ex- 
pected to work well also close to a quantum phase tran- 
sition (T c = 0) driven by varying h. 

III. ID INVERSE-SQUARE FERROMAGNET 

As a non-trivial demonstration of the method, a fer- 
romagnetic chain with interactions decaying as 1/r 2 is 
considered next. The interaction is summed over all i, j 
in l|T]l. i.e., each pair is counted twice. Periodic bound- 
ary conditions are used. includes both distances in 
the periodic system, i.e., 

Ji] = J]i = h (f^tf + (N-\i-j\r) ' (18) 

where J sets the over-all energy scale. 

The classical 1/r 2 Isin g ch ain has been the subject of 
numerous studies [T^, Hj, HE \M& 113 • The long-range in- 
teraction allows for a finite-T phase transition even in one 
dimension. The transition is of an unusual kind, with the 
correlation length exponent v = oo, and a discontinuous 
jump in the magnetization at T c . It can be thought of 
as a one-dimensional analogue of the Kosterlitz-Thouless 
transition, with the topological excitations being kink 
solitons [14|. The model is also important because it 
can be mapped onto the Kondo problem |13| . 

For small h/J, one can expect a behavior similar to 
the classical case, i.e., a finite-T phase transition to a 
ferromagnetic state. For h — ► oo the system becomes 
disordered, and there should therefore be a finite h c for 
which the system undergoes a quantum phase transition 
(i.e., T c — 0). For h < h c , T c > and one can expect 
the same universality class as in the classical case, since 
the quantum fluctuations become irrelevant at T c . Here 
only a single field-strength h/J — 0.5 is considered. The 
simulations show that T c > in this case. 

The model is invariant with respect to flipping all 
spins, which means that for any finite system the average 
magnetization vanishes. The squared magnetization, 

M2 = ((^? < )) ! (19) 

is therefore calculated. Results for M 2 with statistical 
errors in the fifth decimal place can easily be obtained 
for systems with several hundred spins (and there are no 
problems in going to considerably larger systems). For 
small systems the results are in perfect agreement with 
exact diagonalization data. 

A "tempering" scheme, where f3 is considered as an ad- 
ditional discretized dimension of the configuration space 




1.30 1.35 1.40 1.45 1.50 
T/J 

FIG. 5: (a) Magnetization squared vs temperature for system 
sizes N = 16 (dotted curve), 32,64,128,256 and 512 (solid 
curves). The statistical errors are smaller then the width of 
the curves, (b) The same quantity on a more detailed scale in 
the intersection region. The points with barely visible error 
bars are the simulation results. The curves are third-order 
polynomial fits. 

|28j . was also implemented in the simulations. Transi- 
tions satisfying detailed balance are carried out between 
neighboring (3 values. This way, results can be obtained 
on a dense temperature grid with much less effort than 
by several fixed-/3 simulations. A temperature spacing 
AT/ J = 0.01 - 0.02 was used. 

Figure |SJa) shows results for systems with N up to 
512. At high temperatures, M 2 decreases with increas- 
ing N, as expected, and there is a slight increase with N 
at low T. The curves intersect at T/J ss 1.4. A discon- 
tinuous magnetization jump at T c in the thermodynamic 
limit implies that M 2 should become size independent 
at T c for sufficiently large N. A notable difference be- 
tween the finite-size behavior of M(T) seen in Fig. Eland 
the magnetization curves for the classical system is that 
in the latter case the curves do not intersect, but the 
infinite-size value M(T C ) is approached with a logarith- 
mic correction ^l- The reason for the different form of 
the finite-size scaling for h > should be clarified. 
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FigureJSJb) shows in more detail the behavior in the re- 
gion where the curves intersect. The point of intersection 
moves slowly towards higher T as N increases, and larger 
N would be needed to extract T c accurately. Based on 
the data presented here T c / J = 1.42 ± 0.01. This can be 
compared with T c (h = 0) rj 1.53 J for the classical model 
|17| . A reduction of T c is expected on account of quan- 
tum fluctuations for h > 0. The quite small reduction 
for h/J — 0.5 is consistent with the T — > magnetiza- 
tion being only slightly reduced from the classical value 
M(0) = 1. It would clearly also be interesting to study 
the quantum phase transition, but that problem is be- 
yond the scope of this paper. 

The high accuracy of these simulations demonstrate 
that the algorithm indeed is very efficient. The computer 
resources used for this work were quite modest; on the 
order of 200 CPU hours on an SGI Origin2000. The scal- 
ing of the CPU time is close to linear in N for the 1 /r 2 
interaction, for which the interaction sum J^J) converges 
rapidly. Only the local updates discussed in Sec. II B 
were used in these simulations. The cluster updates have 
been tested as well and improve the performance . The 
quantum-cluster update should be particularly useful for 
studying the quantum phase transition, where there will 
be a broad distribution of the sizes of the clusters con- 
structed in this update. 

IV. DISCUSSION 

A new approach to long-range interacting quantum 
models has here been developed within the framework of 
transverse Ising models. It is important to note that the 
technique can also be generalized to other ty pes of sys- 
tems, with the usual caveat of sign problems 24]. What 
is particular about the Ising interaction is that it can be 
written so that a spin-spin term either gives zero or a 
constant when acting on an arbitrary basis state. This 
is what is needed in order to reduce the interactions to 
local constraints in the SSE formalism. However, the 
algorithm can easily be modified to cases where the di- 
agonal interaction can take several non-zero values. The 
first modification is in the diagonal update. For the Ising 
model, the probability of selecting a given bond JSJ) is 
given by a matrix element corresponding to the spin pair 
being in a configuration energetically favored by the in- 



teraction. If the spins are in a non-favored configuration 
(correspinding here to a vanishing matrix element) the 
update is simply rejected. In the general case, the prob- 
ability to use in this update should correspond to the 
largest diagonal matrix element on a given bond, and 
if the actual configuration corresponds to a smaller ma- 
trix clement the update should be accepted only with 
a probability reflecting this smaller value (i.e., the ra- 
tio between the actual value and the largest matrix ele- 
ment). The quantum-cluster update can be modified by 
using ideas developed within the "directed-loop" algo- 
rithm [2i| . For example, there could be 4-particle vertex 
processes where the whole vertex is not necessarily re- 
versed as in Fig. Ufa). The process could instead either 
go straight through the vertex (modifying the vertex only 
at the entrance and exit legs) or "bounce" back without 
modifying the vertex at all. The details of how this is 
done in practice will of course depend on the types of 
diagonal and off-diagonal terms in the Hamiltonian. The 
main point to note is that in the SSE approach all the 
information needed to update the vertices is contained 
in the vertices themselves, which are always local and 
can be generated in the diagonal update based purely on 
local decisions. 

The transverse Ising simulation algorithm has here 
been tested on a one-dimensional model with long-range 
interactions decaying as 1/r 2 . The program requires al- 
most no modifications for higher-dimensional systems, 
and random interactions are also very easy to implement. 
Future studies will have to address how well the method 
works in practice for a variety of systems that are more 
challenging because of frustrated interactions, long-range 
frustrated interactions, or even randomly frustrated long- 
range interactions. For short-range interactions, it would 
also be interesting to see how the SSE quantum-cluster 
method constructed here compares to the transverse Ising 
cluster method previously developed for continuous-time 
worldline simulations [3. 
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